function Output_struc = solve_equi(inputs_struc,fns_struc)

% This function takes all the parameters of the model as given and solves
% for the equilibrium. In particular, it solves for the wage premium which
% clears labor markets and the amount of entry into each quality which
% ensures the zero profit condition holds for each quality level.
% Solving for the free entry conditions takes place in two steps. It is
% possible to first construct a good starting guess for the M_q's by using
% a bisection over 1 dimension combined with a approximation for the ratio
% of M_q's for different q's. Having constructed this approximation, the
% true M_q's are solved for using this as the intitial guess and solving
% for multiplicative deviations around this.

eta = inputs_struc.eta; sigma = inputs_struc.sigma;

if eta == 1/(sigma-1)
    Output_struc = solve_equi_NoLoV(inputs_struc,fns_struc);
else
    Output_struc = solve_equi_LoV(inputs_struc,fns_struc);
end


end


function Output_struc = solve_equi_NoLoV(inputs_struc,fns_struc)

w_U = inputs_struc.w_U; theta_q_i = inputs_struc.theta_q_i;
A_q_i = inputs_struc.A_q_i; sigma_us = inputs_struc.sigma_us; p_z = inputs_struc.p_z;
z_input_q_i  = inputs_struc.z_input_q_i; sigma = inputs_struc.sigma;
N = inputs_struc.N; pdf_A_q_i  = inputs_struc.pdf_A_q_i; eta = inputs_struc.eta;
prod_grid_size=inputs_struc.prod_grid_size; q = inputs_struc.q; a_q = inputs_struc.a_q; 
h = inputs_struc.h;f_q = inputs_struc.f_q; tol = inputs_struc.tol; entry_sh_S_q = inputs_struc.entry_sh_S_q;
no_q = inputs_struc.no_q; delta = inputs_struc.delta; beta = inputs_struc.beta;
A_S = inputs_struc.A_S;


u_fn = fns_struc.u_fn; MC_L_fn = fns_struc.MC_L_fn; L_U_fn = fns_struc.L_U_fn;
L_S_fn = fns_struc.L_S_fn;

% Doing a consistency check on pi. Do I get back the same pi that I put in.
% If not, use the new pi computed and repeat the loop (what ensures its
% convergence?)
iter_pi = 0;
error_pi = 1;
pi_new = 0;

while error_pi>tol && iter_pi<30
    
    iter_pi = iter_pi + 1;
    pi_old = pi_new;

    iter_w = 0;
    error_w = 1;
    min_w = 0;
    max_w = 10;

    while error_w>tol && iter_w<30

        iter_w = iter_w + 1;
        w_S = (min_w + max_w)/2;

        MC_L_q_i = MC_L_fn(w_U,w_S,theta_q_i,A_q_i,sigma_us,A_S);
        p_q_i = (MC_L_q_i + p_z*z_input_q_i)*sigma/(sigma-1);
        P_q_norm = (sum((p_q_i .^(1-sigma)).*pdf_A_q_i ,2)).^(1/(1-sigma));
        P_q = P_q_norm;

        % Share of unskilled and skilled who buy different qualities
        u_q_U = u_fn(q,P_q,w_U,a_q,pi_old);
        rho_q_U = exp(u_q_U)./(sum(exp(u_q_U)));
        u_q_S = u_fn(q,P_q,w_S,a_q,pi_old);
        rho_q_S = exp(u_q_S)./(sum(exp(u_q_S)));

        % Total Demand for a particular quality. 
        D_q = N*(1-h)*rho_q_U.*((w_U+pi_old)./P_q) + N*h*rho_q_S.*((w_S+pi_old)./P_q);

        % Free Entry
        entry_cost_q = (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q);
        M_q_solve = 1/(1-beta*(1-delta))*1/(sigma)*D_q./(entry_cost_q).*P_q;


        d_q_i_temp = M_q_solve.^(eta - sigma/(sigma-1)).*D_q.*(P_q_norm.^sigma);
        d_q_i  = repmat(d_q_i_temp,1,prod_grid_size)./(p_q_i.^sigma);

        % Labor demand for manufacturing production
        l_U_q_i = L_U_fn(d_q_i,A_q_i,theta_q_i,w_U,w_S,sigma_us,A_S);
        l_S_q_i = L_S_fn(l_U_q_i,theta_q_i,w_U,w_S,sigma_us,A_S);

        L_U_q = M_q_solve.*sum(l_U_q_i.*pdf_A_q_i,2);
        L_S_q = M_q_solve.*sum(l_S_q_i.*pdf_A_q_i,2);

        % Labor demand for entry
        n_q = M_q_solve*delta;
        L_entry_U_q = f_q.*(1-entry_sh_S_q).*n_q;
        L_entry_S_q = f_q.*entry_sh_S_q.*n_q;

        % Labor demand from homogenous goods sector
        D_Z_U = 0;%N*(1-h)*rho_B_U*(w_U - p_B)/p_z + N*(1-h)*rho_G_U*(w_U - p_G)/p_z;
        D_Z_S = 0;%N*h*rho_B_S*(w_S - p_B)/p_z + N*h*rho_G_S*(w_S - p_G)/p_z;
        D_Z = 0; %D_Z_U + D_Z_S + D_B*z_input_B + D_G*z_input_G;
        L_Z_U = 0; % L_U_fn(D_Z,A_Z,theta_Z,w_U,w_S,sigma_us,A_Z_U,A_Z_S);
        L_Z_S = 0; %L_S_fn(L_Z_U,theta_Z,w_U,w_S,sigma_us,A_Z_U,A_Z_S);

        % Free Entry
        entry_cost_q = (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q);
        profits_npv = 1/(1-beta*(1-delta))*sum((p_q_i - MC_L_q_i).*d_q_i.*pdf_A_q_i,2);
        free_entry_q = log(profits_npv./entry_cost_q);

        % Labor market clearing
        L_U_clearing = log((sum(L_U_q) + sum(L_entry_U_q))/(N*(1-h)));
        L_S_clearing = log((sum(L_S_q) + sum(L_entry_S_q))/(N*h));
        error_w = abs(L_S_clearing);


        if L_S_clearing>0
            min_w = w_S;
        else
            max_w = w_S;
        end
    end
    
    Pi_q = M_q_solve.*sum((p_q_i - MC_L_q_i).*d_q_i.*pdf_A_q_i,2) - (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q).*n_q;
    Pi_q_test = (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q).*M_q_solve*(1-beta*(1-delta) - delta);
    Pi_q_test2 = max(abs(Pi_q_test./Pi_q - 1));
    pi_new = sum(Pi_q)/N;
    
    error_pi = abs(pi_new*100-pi_old*100);

end

Output_struc.w_S =w_S ;
Output_struc.MC_L_q_i=MC_L_q_i;
Output_struc.rho_q_S=rho_q_S;
Output_struc.d_q_i=d_q_i;
Output_struc.L_U_q=L_U_q;
Output_struc.L_entry_S_q=L_entry_S_q;
Output_struc.L_U_clearing=L_U_clearing;
Output_struc.p_q_i=p_q_i;
Output_struc.P_q=P_q;
Output_struc.P_q_norm=P_q_norm;
Output_struc.D_q=D_q;
Output_struc.l_U_q_i=l_U_q_i;
Output_struc.L_S_q=L_S_q;
Output_struc.free_entry_q=free_entry_q;
Output_struc.L_S_clearing=L_S_clearing;
Output_struc.rho_q_U=rho_q_U;
Output_struc.M_q=M_q_solve;
Output_struc.n_q=n_q;
Output_struc.l_S_q_i=l_S_q_i;
Output_struc.L_entry_U_q=L_entry_U_q;
Output_struc.pi=pi_old;

end



function Output_struc = solve_equi_LoV(inputs_struc,fns_struc)


w_U = inputs_struc.w_U; theta_q_i = inputs_struc.theta_q_i;
A_q_i = inputs_struc.A_q_i; sigma_us = inputs_struc.sigma_us; p_z = inputs_struc.p_z;
z_input_q_i  = inputs_struc.z_input_q_i; sigma = inputs_struc.sigma;
N = inputs_struc.N; pdf_A_q_i  = inputs_struc.pdf_A_q_i; eta = inputs_struc.eta;
prod_grid_size=inputs_struc.prod_grid_size; q = inputs_struc.q; a_q = inputs_struc.a_q; 
h = inputs_struc.h;f_q = inputs_struc.f_q; tol = inputs_struc.tol; entry_sh_S_q = inputs_struc.entry_sh_S_q;
no_q = inputs_struc.no_q; delta = inputs_struc.delta; beta = inputs_struc.beta;
A_S = inputs_struc.A_S;

u_fn = fns_struc.u_fn; MC_L_fn = fns_struc.MC_L_fn; L_U_fn = fns_struc.L_U_fn;
L_S_fn = fns_struc.L_S_fn;

% Doing a consistency check on pi. Do I get back the same pi that I put in.
% If not, use the new pi computed and repeat the loop (what ensures its
% convergence?)
iter_pi = 0;
error_pi = 1;
pi_new = 0;

while error_pi>tol && iter_pi<30
    
    iter_pi = iter_pi + 1;
    pi_old = pi_new;

    iter_w = 0;
    error_w = 1;
    min_w = 0;
    max_w = 10;

    while error_w>tol && iter_w<30

        iter_w = iter_w + 1;
        w_S = (min_w + max_w)/2;

        MC_L_q_i = MC_L_fn(w_U,w_S,theta_q_i,A_q_i,sigma_us,A_S);
        p_q_i = (MC_L_q_i + p_z*z_input_q_i)*sigma/(sigma-1);
        P_q_norm = (sum((p_q_i .^(1-sigma)).*pdf_A_q_i ,2)).^(1/(1-sigma));


        % For a given w_S, we now try to solve the free entry conditions i.e.
        % the M_q's which solve the free entry conditions. First construct an
        % approximate solution using 1 dimensional bisection and ratios of
        % M_q's. Then use this as a starting point to fsolve for M_q's.
        % The next 4 lines get used in constructing different M_q's relative to
        % the base q's' M_q. The base q's M_q is found by bisection.
        g = sum(exp(a_q).*((w_U./P_q_norm).^q))./sum(exp(a_q).*((w_S./P_q_norm).^q));
        M_power = sigma/(sigma-1) - eta + (q+1)*(eta - 1/(sigma-1));
        M_ratio_1 = (f_q./f_q(1,1)).*(exp(a_q(1,1))./exp(a_q)).*(P_q_norm(1,1).^(-q(1,1)))./(P_q_norm.^(-q));
        M_ratio_2 = ((1-h).*((w_U+pi_old).^(1+q(1,1))) + h*((w_S+pi_old).^(1+q(1,1))).*((w_U+pi_old).^(mean(q)))./((w_S+pi_old).^(mean(q))))./...
            ((1-h).*((w_U+pi_old).^(1+q)) + h*((w_S+pi_old).^(1+q)).*((w_U+pi_old).^(mean(q)))./((w_S+pi_old).^(mean(q))));

        % Find an approximate solution for free entry co
        iter_M_q = 0;
        error_M_q = 1;
        min_M_q = 0;
        max_M_q = N;

        while error_M_q(1,1)>tol && iter_M_q<30

            iter_M_q = iter_M_q + 1;
            M_q_guess = (min_M_q + max_M_q)/2;
            % Approximation for different M_q's for a given base M_q. This is
            % based on the equation labeled "entry ratios" in the lyx file
            % "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx".
            M_q_approx = ((M_q_guess.^(M_power(1,1)))./M_ratio_1./M_ratio_2).^(1./M_power);

            P_q = P_q_norm.*(M_q_approx.^(eta - 1/(sigma-1)));

            % Share of unskilled and skilled who buy different qualities
            u_q_U = u_fn(q,P_q,w_U,a_q,pi_old);
            rho_q_U = exp(u_q_U)./(sum(exp(u_q_U)));
            u_q_S = u_fn(q,P_q,w_S,a_q,pi_old);
            rho_q_S = exp(u_q_S)./(sum(exp(u_q_S)));

            % Total Demand for a particular quality. 
            D_q = N*(1-h)*rho_q_U.*((w_U+pi_old)./P_q) + N*h*rho_q_S.*((w_S+pi_old)./P_q);

            % Demand for each point on grid.
            d_q_i_temp = M_q_approx.^(eta - sigma/(sigma-1)).*D_q.*(P_q_norm.^sigma);
            d_q_i  = repmat(d_q_i_temp,1,prod_grid_size)./(p_q_i.^sigma);

            % Free Entry
            entry_cost_q = (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q);
            profits_npv = 1/(1-beta*(1-delta))*sum((p_q_i - MC_L_q_i).*d_q_i.*pdf_A_q_i,2);
            free_entry = log(profits_npv./entry_cost_q);

            error_M_q = abs(free_entry);

            if free_entry(1,1)<0
               max_M_q = M_q_guess;
            else
                min_M_q = M_q_guess;
            end


        end


        % Create a new structure which only has inputs which are needed for
        % solving the free entry condition exactly. I could have used the
        % inputs structure itself but not changing things here due to history
        % dependence!
        inputs_Mq_struc.eta               = eta;               
        inputs_Mq_struc.sigma             = sigma;
        inputs_Mq_struc.w_U               = w_U;             
        inputs_Mq_struc.w_S               = w_S;        
        inputs_Mq_struc.N                 = N;       
        inputs_Mq_struc.h                 = h;                 
        inputs_Mq_struc.entry_sh_S_q        = entry_sh_S_q;
        inputs_Mq_struc.P_q_norm          = P_q_norm;          
        inputs_Mq_struc.q                 = q;                 
        inputs_Mq_struc.a_q               = a_q;               
        inputs_Mq_struc.f_q               = f_q;               
        inputs_Mq_struc.pdf_A_q_i         = pdf_A_q_i;         
        inputs_Mq_struc.u_fn              = u_fn; 
        inputs_Mq_struc.prod_grid_size    = prod_grid_size;
        inputs_Mq_struc.p_q_i    = p_q_i;
        inputs_Mq_struc.MC_L_q_i    = MC_L_q_i;
        inputs_Mq_struc.M_q_approx = M_q_approx; 
        inputs_Mq_struc.delta    = delta;
        inputs_Mq_struc.beta = beta;
        inputs_Mq_struc.pi_old = pi_old;


        options_m = optimset('Display','off');
        % Solve for free entry exactly. Note that the fsolve does not take the
        % M_q's as the 'x' variable to put to 0. Rather it solves for the
        % multiplicative deviations of the M_q's from the starting guess. This
        % is because the M_q's have very different scales which leads to
        % inaccuracies in the gradient based solvers. The good starting guess
        % constructed earlier means that the multiplicative deviations lie
        % around 1 which can be solved for very accurately.
        scale_Mq_solve = fsolve('M_q_zero',ones(no_q,1),options_m,inputs_Mq_struc);
        free_entry = M_q_zero(scale_Mq_solve,inputs_Mq_struc);
        % Using the solved for multiplicative deviations to figure out entry.
        M_q_solve = M_q_approx.*scale_Mq_solve;

        P_q = P_q_norm.*(M_q_solve.^(eta - 1/(sigma-1)));

        % Share of unskilled and skilled who buy different qualities
        u_q_U = u_fn(q,P_q,w_U,a_q,pi_old);
        rho_q_U = exp(u_q_U)./(sum(exp(u_q_U)));
        u_q_S = u_fn(q,P_q,w_S,a_q,pi_old);
        rho_q_S = exp(u_q_S)./(sum(exp(u_q_S)));

        % Total Demand for a particular quality. 
        D_q = N*(1-h)*rho_q_U.*((w_U+pi_old)./P_q) + N*h*rho_q_S.*((w_S+pi_old)./P_q);


        d_q_i_temp = M_q_solve.^(eta - sigma/(sigma-1)).*D_q.*(P_q_norm.^sigma);
        d_q_i  = repmat(d_q_i_temp,1,prod_grid_size)./(p_q_i.^sigma);

        % Labor demand for manufacturing production
        l_U_q_i = L_U_fn(d_q_i,A_q_i,theta_q_i,w_U,w_S,sigma_us,A_S);
        l_S_q_i = L_S_fn(l_U_q_i,theta_q_i,w_U,w_S,sigma_us,A_S);

        L_U_q = M_q_solve.*sum(l_U_q_i.*pdf_A_q_i,2);
        L_S_q = M_q_solve.*sum(l_S_q_i.*pdf_A_q_i,2);

        % Labor demand for entry
        n_q = M_q_solve*delta;
        L_entry_U_q = f_q.*(1-entry_sh_S_q).*n_q;
        L_entry_S_q = f_q.*entry_sh_S_q.*n_q;

        % Labor demand from homogenous goods sector
        D_Z_U = 0;%N*(1-h)*rho_B_U*(w_U - p_B)/p_z + N*(1-h)*rho_G_U*(w_U - p_G)/p_z;
        D_Z_S = 0;%N*h*rho_B_S*(w_S - p_B)/p_z + N*h*rho_G_S*(w_S - p_G)/p_z;
        D_Z = 0; %D_Z_U + D_Z_S + D_B*z_input_B + D_G*z_input_G;
        L_Z_U = 0; % L_U_fn(D_Z,A_Z,theta_Z,w_U,w_S,sigma_us,A_Z_U,A_Z_S);
        L_Z_S = 0; %L_S_fn(L_Z_U,theta_Z,w_U,w_S,sigma_us,A_Z_U,A_Z_S);

        % Free Entry
        entry_cost_q = (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q);
        profits_npv = 1/(1-beta*(1-delta))*sum((p_q_i - MC_L_q_i).*d_q_i.*pdf_A_q_i,2);
        free_entry_q = log(profits_npv./entry_cost_q);

        % Labor market clearing
        L_U_clearing = log((sum(L_U_q) + sum(L_entry_U_q))/(N*(1-h)));
        L_S_clearing = log((sum(L_S_q) + sum(L_entry_S_q))/(N*h));
        error_w = abs(L_S_clearing);


        if L_S_clearing>0
            min_w = w_S;
        else
            max_w = w_S;
        end
    end

    Pi_q = M_q_solve.*sum((p_q_i - MC_L_q_i).*d_q_i.*pdf_A_q_i,2) - (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q).*n_q;
    Pi_q_test = (f_q*w_U.*(1-entry_sh_S_q) + f_q*w_S.*entry_sh_S_q).*M_q_solve*(1-beta*(1-delta) - delta);
    Pi_q_test2 = max(abs(Pi_q_test./Pi_q - 1));
    pi_new = sum(Pi_q)/N;
    
    error_pi = abs(pi_new*100-pi_old*100);
    
end


Output_struc.w_S =w_S ;
Output_struc.MC_L_q_i=MC_L_q_i;
Output_struc.rho_q_S=rho_q_S;
Output_struc.d_q_i=d_q_i;
Output_struc.L_U_q=L_U_q;
Output_struc.L_entry_S_q=L_entry_S_q;
Output_struc.L_U_clearing=L_U_clearing;
Output_struc.p_q_i=p_q_i;
Output_struc.P_q=P_q;
Output_struc.P_q_norm=P_q_norm;
Output_struc.D_q=D_q;
Output_struc.l_U_q_i=l_U_q_i;
Output_struc.L_S_q=L_S_q;
Output_struc.free_entry_q=free_entry_q;
Output_struc.L_S_clearing=L_S_clearing;
Output_struc.rho_q_U=rho_q_U;
Output_struc.M_q=M_q_solve;
Output_struc.n_q=n_q;
Output_struc.l_S_q_i=l_S_q_i;
Output_struc.L_entry_U_q=L_entry_U_q;
Output_struc.pi=pi_old;

end